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^ I Abstract. We have developed a new numerical technique for simulat- 

ing dusty-gas flows. Our unique code incorporates gas hydrodynamics, 
self-gravity and dust drag to follow the dynamical evolution of a dusty-gas 
Q> ■ medium. We have incorporated several descriptions for the drag between 

i gas and dust phases and can model flows with submillimetre, centimetre 

^ ■ and metre size "dust". We present calculations run on the APAC^ super- 

computer following the evolution of the dust distribution in the pre-solar 
nebula. 
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^1 1. Introduction 

Up until recently we had only one observation to test our theories of planet 
^ . formation against - our own Solar System. Now however with new planets and 

^ I solar systems being both identified and parameterised at a rate approaching 

one a month, the observational constraints are much tighter and our lack of 
understanding of many aspects of the planet formation process is all too obvious. 
At the most basic level we know that micron size grains of dust in the pre-solar 
nebula clump and coagulate together to form planets, objects 10^^-10^^ times 
larger. Planet formation is a multi-stage process, taking us from dust grain to 
boulder to planetesimal to planetary embryo. Analytical arguments (Goldreich 
& Ward 1973) have presented us with constraints on the time scales for each 
stage but little more. It is the very first stage of the process that we are concerned 
with in this paper - from micron scale dust to metre sized boulders. 

Theoretical models have changed much in recent years and the simple pic- 
ture of a thin dust layer accumulating at the disk midplane, becoming gravita- 
tionally unstable, and breaking into planetesimals (Safronov 1969; Goldreich & 
Ward 1973) now seems unlikely. Recently Goodman & Pindor (2000) showed 
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that turbulent drag causes radial instabilities in the dust layer, even if the disk 
self-gravity is negligible. In their steady state models, grains in a uniform dust 
layer experience radial drift as expected. However, for perturbed disks they pre- 
dict that over-dense rings form within an orbital period, with the ring thickness 
similar to the thickness of the dust layer. These rings eventually collapse into 
planetesimals in the kilometre size range. 

In this paper we present the first three-dimensional numerical simulations 
that include the effects of hydrodynamical forces, self-gravity and gas drag upon 
an evolving dusty gas disk. We describe a new numerical code, based upon the 
smoothed particle hydrodynamics (SPH) technique which uses a collection of 
particles to approximate a fluid. At present we run simulations with uniform 
grain size and do not allow dust particles to have individual or time varying 
chemical and physical properties. Adding this level of sophistication to the 
model would require only minor changes to the code. 



2. The Planet Building Code 

We have developed a new parallel two-phase (dust & gas) hydro -t- tree code by 
merging a dusty gas SPH code (Maddison 1998) with an MPI parallel Hashed 
Oct Tree N-body -|- SPH code (Humble 1999). Following Monaghan Sz Kocharyan 
(1995), the dusty gas is approximated by two inter-penetrating flows that inter- 
act via a drag force. 

The particles are tagged as being gas or dust and are then allowed to interact 
according to gas-gas, gas-dust and dust-dust interactions. Only gas particles 
feel a pressure force and are affected by viscosity. The gas-dust combination feel 
the drag force and also the mixed term seen in both the dust and gas particle 
acceleration equations. The only dust-dust interaction is through gravity. 

The two fluids are coupled by gravity and drag, and the equation of motion 
is given by: 

^ = -i VP+-(v,-v,)-V$, 

dt Pg Pg 

!^ = _ivP-f(v,-v,)-V.f, 

dt Pd Pd 

where the terms on the right hand side are gas pressure, drag force and gravity 
respectively. The density p is the mass density per unit volume, while p is the 
local density (related via p = Qp, where Q is the volumetric void fraction). The 
functional form of K depends upon the drag regime (e.g. Epstein or Stokes) 
being considered. In the simulations we present, we have implemented Epstein 
drag which is appropriate for protoplanetary disks. Thus K = (/90(7dragc)/?^dusti 
where c the local sound speed, Cdrag the drag coefficient, and rdust the dust grain 
size. 

The drag term is calculated using a pairwise implicit backward Euler scheme, 
while the hydrodynamics and self-gravity are solved explicitly. Operator split- 
ting is used to integrate particle orbits under the influence of the un-softened 
gravity of the central star. The code has been tested for stability and accuracy 
with standard periodic box simulations to ensure that the drag terms are correct. 
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Figure 1. The 4 left hand panels show final gas (dots) and dust 
(contours) density distribution for 12. 5K gas + 12. 5K dust models of 
various grain sizes. Five contour levels (10^^^ to 10^^'' g/cm"') arc 
plotted. The 2 right hand panels show the final dust surface density 
distributions of all 6 models (initial distribution labelled to). 



and collapsing sphere simulations to test the gravity. The code scales linearly to 
more than 32 processors. For simplicity we assume that the dust grains do not 
evaporate or coagulate and that the gas does not condense. We take the dust 
grains to be incompressible. 



3. The Simulations 



We are particularly interested in the effects of grain size upon disk morphology 
and global disk dynamics. We present the results of 6 simulations of gas disks 
laden with dust of different grain sizes. All simulations start with a prolate 
spheroid of rotating, self-gravitating gas (with a star at the centre) that spins 
down to near-equilibrium. The dust is then added (overlaid on the 3D flared gas 
disk) and the simulations are then followed for approximately 10^ years (which 
corresponds to 11 orbits at 100 AU). 

The model parameters used were: = 1.0 Mq, M^isk = 0.01 M^,, Mjust = 
O.OlMdisk and -Rdisk = 100 AU. The grain parameters were: pd = 2.4g/cm'^ and 
Td = 1, 10~^, 10~^, 10~^^ 10~^ and 10~^m. The a-disk parameters used were: 
c{R) = co{R/ 100 AU)-^/^, T{R) oc {R/ 100 AU)-^/^ , H/R = 0.1 at 100 AU, and 
7 = 5/3 with an isothermal equation of state. Low resolution runs used 2 x 12500 
particles, and high resolution runs used 2 x 125000 particles. The self-gravity 
softening e = 1.0 AU, and the SPH artificial viscosity parameters cksph = 0.1 
and /3sPH = 0.0. 
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4. Discussion 

For large (10m) and small (micron) dust sizes, we expect that the dust distribu- 
tion will stay close to the initial flared disk. The largest grains (not shown) are 
weakly coupled to the gas, and if started in Keplerian motion, they will remain 
there. On the other hand, the tiny grains are so strongly coupled to the gas 
that they are essentially co-moving (on the timescales we examine). For both 
extremes, we see little evolution of the dust distribution. It is for the regime 
0.1 mm < rdust < 1 m that the most significant disk evolution occurs. 

In the r — z plots of figure 1, significant deviation from the initially flared 
disk occurs in the inner regions of the Im and 10cm plots, in the mid regions 
(from r = 0.6 to r = 0.9) of the 1cm plot, and in the outer regions of the 1mm 
and smaller plots. In these regions, the dust is moving at close to Keplerian 
speeds, whereas the gas is (as always) sub-Keplerian. As the velocity difference 
is sustained, the Epstein drag is optimal, and the energy loss rate is high — 
allowing a thin layer of dust to form in the midplane and migrate radially. 
These thin dense dust disks are those that Goodman &; Pindor suggest have 
global turbulent instability modes. 

While the 10cm and 1cm dust exhibits the highest surface density (top right 
of figure 1), the volumetric density is largest in the inner regions of the Im and 
10cm disks. Therefore these size ranges are probably the most interesting from 
a planet formation viewpoint. 

The Lagrangian nature of the code means that it is trivial to add empirical 
grain growth models, and to follow the grain temperature and density histories 
and hence generate chemical compositions. The equations of state and drag 
term can easily be altered on a per-region or per-particle basis to account for 
local disk conditions. 
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